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A nonuniform sheet of ions generated at time t = diffuses in a cylindrical ion diffusion tube con- 
taining a nonreacting neutral species flowing with parabolic velocity distribution. Calculation of the 
on-axis ion density at a point z downstream as a function of time t is reduced to a single numerical inte- 
gration for each (z, t) involving some functions which have been computed once for all. An example 
is given showing the effect of the velocity distribution compared with a uniform flow with the same 
flow rate. The results appear to be corroborated by experiment. 
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1 . Introduction 

The analysis presented here is currently being used to predict the behavior of a diffusing ion 
sheet in a flow tube. In the experiment [5], 1 an inert buffer gas, usually helium, is injected into a 
long flow tube about 8 cm in diameter. Distances along the tube are measured from a point z = 
downstream at which the flow has become essentially laminar, with a mean velocity of 8(10) 3 cm/s 
and a pressure of 0.4 torr. A thin sheet of ions is generated at z = Q at time t= by a pulsed beam of 
electrons incident perpendicular to the axis of the tube. The pulse width is variable, but 100 jjls 
might be considered characteristic. The relatively high pressure in the tube insures that the ions 
created will be thermalized at gas temperature in a few centimeters. The ion sheet is carried down- 
stream by the flowing gas with a pulse shape determined by diffusion and the nonuniform velocity 
profile. 

The sampling port of a mass spectrometer is located on the axis at a variable distance z down- 
stream from the electron gun, where z is of the order of 100 cm. The detector is gated with a gate 
width of 100 jus. Thus the ion arrival spectrum as a function of delay time between the termination 
of the excitation pulse and the initiation of the detection pulse is measured. Using experimentally 
measured values for the diffusion coefficient Z), the computed arrival spectrum agrees well with the 
actual measured spectrum, indicating no serious errors are present in the formulation and solution 
of the transport equation. 

Letting (r, 6) be coordinates in the cross section of the tube, the problem is to determine the 
time-dependent ion density G(r, 6, z, t) on the axis (r=0) at a given distance z downstream 
from the point at which the ions were generated as a thin sheet with given (r, 6) -dependence. 
The background flow is assumed to be laminar and steady, with a parabolic velocity distribu- 
tion, containing perhaps a neutral reactant species of uniform density F everywhere. With a 
reaction rate k, G obeys the differential equation: 
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where o§? is the Laplacian in cylindrical coordinates. Here v(r) is the parabolic velocity distribution: 

v(r) = 2v(l-($y (2) 

where V is the mean velocity and a is the tube radius. At the wall of the tube, the ion density is 
assumed to be negligible: 



G(a, 0,2, t) = 0. 
We shall take for the given initial ion density distribution: 

G(r, 0,z, 0) = P(r, 0)8(z). 



(3) 



(4) 



We notice that if kF = 0, the resulting solution G (r, 0, 2, t) gives the general solution G for 
kF t^ as follows: 

G(r, 0, 2, *) = G (r, 0, 2, ^" KFr , (5) 

so that we may take kF = in (1) without loss of generality. 

2. Solution by Integral Transform 

We shall expand G by Laplace transform in t, Fourier transform in 2, and an azimuthal expan- 
sion in 0. 

C(r, 0, 2, t) = -^-- ? f X rfAre** r^ f °° dse st Y e' n "F w (r; A, 5), 
(2t7) 2 J- m 2ttj Jc-ix „~ x 



Too foo T27T 

F„(r; A, s)= &e-** fifar*' dOe- in0 G(r, 0, 2, *). 

J-x Jo Jo 



(6) 



Then 



^hei 



d 2 F n 1 SF n 

dr 2 r dr 



n 2 72 s + ikv(r) 



Z) 



Fn = -pPn(r), 



■r 
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We introduce the variables 



P„(r) = rf0P(r,0). 

Jo 

F H (r)=T7»e-*' 2 #„(Tj), 
P»(r) = 9„(tj), 
T7 = 7r 2 , 



(8) 



(9) 



where 
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f-^t-a)- 



(10) 



We obtain the inhomogeneous confluent hypergeometric equation [la] 

1 qn{y])e^ 2 



V m+ (l + 2n-7i)H' H -a H H„-- 



vhere 



1 1 

a n = n + - + — * 

2 4y 



k 2 + 



- 4yZ) 77" 



s + 2iA:H 



D 



The homogeneous equation has the solutions [lb, lc] 

4>(a„, 1 + 2n; 17), i//(a„, 1 + 2/z; 17), 

with Wronskian 

e^r(l + 27i) 



w- 



r? 1+2n r(a„) 



(11) 



(12) 



(13) 



(14) 



As must be expected, only the lowest azimuthal modes, n = 0, contribute to the axial value 
r=0 of the ion density. The solution satisfying smoothness at 77 = and vanishing at r) = ya 2 is, 
for n = Q: 



Fo(r) 



QSo) 
4yZ) 



jgQyg 2 ) 



(/)(a , 1; ya 



77 c/)(a , 1 



;T7)-^(r7)}, 



where 



&(ri) = 6(aQ, 1; 77) d«e 2 qo(t)ib(ao, 1; r) — iMao, 1; 17) d«e 2 , 

Jo Jo 

Onr=0, 



qo(t)(f)(ao, 1; f). 



If we write 



/ 7 o(0)=^r(a„)^(ya 2 )/</.(5„, 1; ya 2 ). 



b= a 2 y 



and 



^(S;oFo(0;i,5)*x(«*) 

for the inverse Laplace transform, (17) gives 

a 2 f -— f A 

x(/;£)=— I (fee 2 9o(6z)=Sf(l! o jr(a ) • t//(ao, 1; bz) — T(a ) 



ift(ao, 1; b)<f>(ao, 1; 6z) 
0(a o , 1; 6) 



We put 



a 2 = 4D 

46D 5 ' T a 2 *' 

v 4 A + w ik, 



(15) 



(16) 



(17) 



(18) 



(19) 



(20) 
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and (19) becomes by translation 

X (t: k) = be-( p +l) T J* dzeVqtfz)££ M • ^(A)^A, 1; bz)-Y(A)^(A, 1; b) T^^f}' (21) 

The J?f _1 of the first term in the bracket can be found in tables [2a], when </> is replaced by the 
proper Wittaker function, or may be found by the residue theorem from the poles of T{A) at 
A = 0, — 1, — 2, . . .. The second term may be evaluated in two parts. First, the poles of T(A) are 
shown in the appendix 1 to give a term which exactly cancels the contribution from the first term 
in the bracket. Second, the zeroes of 4>{A, 1; b), which we shall designate as a m (b): (m =1,2,. . .), 

<j>(a m (b\ 1; 6) = 0, (22) 

give by the residue theorem 

x{t; k ) = e-^T(b, r), r(6, r)= J T m (b, r), 

r,„(6, r) — bfl ^""'^ *+ ■ J *ggs« (23) 
0,a: o (a w , 1; b) ££ (X.!)- 

The last factor in T m {b, r) comes from integration with respect to z using the power series for 
</>(a w , 1; bz): 



*E(b;Po)= f dze~\ z z K P»(aVz) 
Jo 



J?K(b; P ) = dze 2 2 z K Po(aVz). (24) 



It is important to note that when we let b — » in the above expression, we find 

f ( "7F# «) "* 2 f ' dzzJo(£ m z)Po(az), 

using the limiting forms for $ in terms of 7 as b — > 0, a m 6 finite: [Id] 



(a, 1; W~ /o(2V6S). 
|6a| < °° 



(25) 



This also gives the location of zeroes for small b: 

a m (b)-*-gJ4b, (26) 

where Jo(^m) = 0, (m= 1, 2, 3, . . .). Indeed, as 6^0, the entire expression (23) becomes 

independent of b, and the inverse Fourier transform involves only the factor e~ PT , from which we 
recover the solution corresponding to a uniform velocity v(r) = 2V. ForP (0 = 27r,this is 



G(o,z, «)=: 


1 


401 2 f«/,(f„)' 


V^Di £ 


v(r) m 2V. 







(27) 
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For z > 0, the { — ik) 112 factor in b is to be interpreted as lim ( — / — iky 12 , (/ > 0), so that the 
path of integration in the /r-plane for the inverse Fourier transform passes below the branch at k = 0: 

- 7T ^ arg k ^ 0. (28) 

We replace k by b through: 

k= Sv b2et *' (29) 

and take the path of integration C in the 6-plane to run from o°e _i T to the origin, Ci, and out to 

ooe"' ?, C 2 . 

Then with 

7 = 2D' ^ = 2^F Z ' (30) 

we have 

G(0, z, /) =^ | <#i6r(6, x)e-" 2 f - I) + ^. 

Putting 6 = ue~ j T on Ci, 4= ae~' 4 on C2, and noting that G is real, we obtain 

G(0, z, = « ^ ea ^ I duue~ J^y u4t+iu ( ~ 4) • r(ae _/ ?, r). 

Z7ray Jo 



(31) 



(32) 



In a particular case of interest, (8y) 2 = 3.7(10) :J , ^ — 0.3, 



T 



0.1. For such values 



4 

expansion methods are not satisfactory, especially in view of the difficulty of getting any simple 
expressions for the zero-curve a m (b) and F(b, r). Numerical integration then seems to be necessary. 



It is convenient to write 



3. Numerical Integration 



C(0, z, = £ Cm(z, t), (33) 



where 



G m {z, t) = Real dime (*y) 2 " "' ( 4 ) ' r m (ue-' 1, r), 

liray Jo 



(34) 



27T(iy 

where the subscript m refers to the terms due the zeroes of (f)(a m (b), 1; ue~ i 4) which, as u —> 
are given in (26). For computational purposes, it is convenient to introduce the following functions: 

JF(u 9 B) =<t>(Blu, I; ac-U), 

W,n(i£, £) =</>, * (fi/M, 1; ue-'f). (35) 
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Treating u as a parameter, the complex zero curve for each m, 

tF(u,B m (u))=0 9 
[a m (b) = B m (u)/u, 
is found by Newton-Raphson iteration 

F(u,BW) 






F, B (u,B®Y 



B\V(0)=--Zle>i. 



(36) 



(37) 



using the expansions [16] 



F(u,B) = ^Q K {u,B). 



(38) 



where 



Qk(u, B) 



B(B + u) . . . (B+(K-l)u) .km 



(K\y 



F, H (u,B)=^Q K (u,B)-[^+-^- u 



+ . . . + 



1 



B+(K-l)u] 



(39) 



with Kf chosen to control relative error in Fh at cutoff. Curves B m (u) for m= 1, 2, 3 were computed 
for ^ u ^ 17 for a set of Gaussian quadrature abscissae on each unit interval in u. 

For fixed u, the number N of zeroes of (/> I , 1; ue~ j 4 I inside a circle radius p centered 

on Bo is easily shown to be 



N 



-$/>■ 



F, B (u,flo-pe'") 
F(u, Bo-pe'"*) ' 



(40) 



Using Gaussian quadrature, it was easily established that no zeroes occur for |B| ^ 10, 
^ u ^ 17, other than those continued from B,„(0). 

The remaining functions necessary to compute T(b, r) are found from the following: 



^(ii) s-rW^a.,, 1; tie **) = % Qk(u 9 B m ) [*(— +K) -2*(1+K) 



(41) 



/„(u; Po)= J <?*(", B m )S K (oe J ; P»), 



(42) 



where i//(z) is the Digamma function [3a] and -A' is denned in (27). t//(z) is computed by asymptotic 
series for |z| > 4, and for \z\ <4 the recursion relation is used to get z outside this circle. 

Clearly the radial distribution Po{r) occurs only in l m (u) via «/#, and does not affect the 
zero-curve location. Then 



r,„(we 4 ,T)=e 4 e - - — 

F,«(ii, B m (u)) 



(43) 
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and the numerical integration in (32) is carried out by Gaussian quadrature on the unit intervals 
of u. 

From the discussion following (24), the distribution for uniform velocity profile v(r) = 2V is 



[G(0, Z, t) J uniform - 



1 



27ray 



Real 



• I duu 

Jo 



e (8-y) 



7^b^ + - 2 (^-?) 



r(0; r). 



4. Results: P (r) = 2tt 

Detailed calculations have been made only for the uniform radial distribution Po(r) = 277 for 
which (24) is readily evaluated by recursion. (In particular, the results apply to a uniform distribu- 
tion P(r, 6) = 1 in (8)). A numerical evaluation for more general Po(r) is readily substituted. For 
typical numbers (r, j?, y) of interest, the first term, m= 1 in (33), dominated the second, m = 2, by 
the factor 10 4 . 

Figure 1 shows B m (u) , the zero-curve, for m= 1, 2. Figure 2 shows 6(0, z, t) for a = 4 cm, 
Z)= 1800 cm 2 /s, z = 42 cm, and a mean velocity F=6800 cm/s for both the parabolic and uniform 
velocity profiles. 

Figure 3, taken from a forthcoming publication by Fehsenfeld [5], shows a comparison of 
the experimental profile (solid points) and the calculated profile (open points) for the values a = 4 cm, 
Z) = 2400 cm 2 /s, z=113 cm and F=6800 cm/s. In this figure, both sets of data are normalized to 
peak values of unity, and the time of the experimental peak is shifted by 0.3 ms. 

Computation of the zero-curve for each m took three minutes, while the curves in figure 1 
took two minutes each on the CDC 3600 with complex arithmetic hardware. A 50-point Gaussian 
integration was used throughout [4]. Fortran programs are available from the author. 



The author wishes to express his thanks to Fred C. Fehsenheld and Arthur L. Schmeltekopf 
for suggesting the problem and providing experimental verification of the results. 
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Figure 1. Zero curves B m (u) for m = 1, 2 (eq 37). 



Figure 2. Comparison of sample ion densities for uni- 
form and parabolic velocity profiles. 
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FIGURE 3. Comparison of normalized experimental and 
calculated profiles with juxtaposed peaks. 
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6. Appendix 

We shall show that contributions to the inverse Laplace transform due to the poles of T(A) 
in the two terms in (21) cancel out. 

From well-known relations [le, 2a], the first term gives: 



bz 



^ A \ bT) T(A)x\>(A, 1; bz)=J?~^ (bz)-Wr(A)W l/2 -A,o(bz) 



e bT (e bT -\)- x exp 



V» T -i/' 



(Al) 



according to tables [2a]. (Note that in the reference [2a], second column, (2 — K — u)t should 
read (2+K— u)t to agree with reference [2b] and the straightforward residue evaluation.) 
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With regard to the second term, let us put 

and assume that this factor does not contribute to Jzf~ l . (Its contributions are treated separately 
in the text.) Then for c > 1 (the nonlogarithmic case) [lc] 

JTft^UWA, 1; b)*(A) =Q(t; b, c) 

sin 7rc [ r(c)r(i+^— c) 

_ bl _ c ^ <b(l + A^2-c; h ) MA) ^ (A3) 

by the definition of <//(^, c; 6) in terms of two independent solutions (f>(A, c; b) and 6 1_c $(l H-/4 — c; 
2 — c; 6) for noninteger c. The second term on the right of (A3) vanishes due to the regularity of 
the integrand, so that by the residue theorem, the poles of V(A) give: 

Vl ' ' ' sinTrc^^ ; m\r(c)T(l-m-c) l '' 1A4j 

m=0 v / v / 

Taking the limit c — > 1 + , 

<f)(— m, c; 6) — * 0(— m, 1; 6), 

r(c) - i, 

(-)"'7T 

sin 7rcr(l — m — c) — > : — , (c—> 1 + ). 

We find then, that c -> 1 + , 

(&; 6, 1)= J <f>(-m, 1; 6)e- m *<D(-ro). (A5) 

m=0 

Replacing 3>(—ra) by (A2), 

(?(*;6,1) = J <M-m,l;6z)e-«*. (A6; 

m = 

By the formula for the generating function for the Laguerre polynomials [3b], 

J <M- >*, 1; &)* m = (1-*) - 1 exp (—); (A7) 

putting jc = e _6r , we obtain the equivalence of the expressions (Al) and (A6). 

(Paper 72B2-267) 
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